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ABSTRACT 

This  paper  describes  a  package  for  analytical  ray  tracing  of 
relatively  simple  optical  systems.  AESOP  (An  Extensible 
Symbolic  Optics  Package)  enables  analysis  of  the  effects  of 
small  optical  element  misalignments  or  other  perturbations. 
(It  is  possible  to  include  two  or  more  simultaneous 
independent  perturbations.)  Wavefront  aberrations  and 
optical  path  variations  can  be  studied  as  functions  of  the 
perturbation  parameters.  The  power  of  this  approach  lies  in 
the  fact  that  the  results  can  be  manipulated  algebraically, 
allowing  determination  of  misalignment  tolerances  as  well  as 
developing  physical  intuition,  especially  in  the  picometer 
regime  of  optical  path  length  variations. 


a 


1.  Introduction 


Modem  high-precision  optical  systems,  such  as  space  astrometric  interferometers  (e.g.  EAME: 
Johnston  et  al.  1997;  POINTS:  Reasenberg  et  al.  1996;  GAIA:  Eoiseau  and  Malbet  1996,  Eoiseau 
and  Shaklan  1996;  Eindegren  and  Perryman  1996;  SIM:  SIM97),  can  require  optical  path 

(-9)  (-12) 

tolerances  in  the  sub-nanometer  {\nm=\Q  m)  to  picometer  (1  pm  =  10  m)  regimes  over 
total  path  lengths  on  the  order  10  m.  Such  tolerances  place  extreme  requirements  on  optical 
analysis  programs.  Two  questions  are  of  paramount  importance:  1)  to  which  specific 
perturbations  is  a  system  most  sensitive?  and  2)  are  there  couplings  between  different 
perturbations  that  produce  high  sensitivities  (i.e.,  are  there  strong  correlations  between 
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perturbation  parameters)?  AESOP  can  be  used  to  answer  these  questions  (Munson,  1993),  as  well 
as  to  develop  physical  intuition  in  the  picometer  optical  path  difference  (OPD)  regime. 


A  common  optical  subsystem  employed  in  astronomical  interferometers  is  a  beam  compressor, 
used  to  convert  a  large  aperture  input  beam  (starlight)  to  a  narrow  output  beam  (~1  cm)  suitable 
for  combining  with  another  such  beam  to  produce  interference  fringes  of  sufficient  visibility.  A 
typical  beam  compressor  consists  of  a  pair  of  confocal  paraboloidal  mirrors,  as  sketched  in  Figure 
1 .  If  perfectly  aligned,  a  flat  input  wavefront  results  in  a  radially  compressed  flat  output 
wavefront.  Misalignment  analysis  of  even  such  a  simple  system  as  this  generally  requires 
resorting  to  numerical  programs.  Usually,  such  programs  are  ill-suited  for  studies  involving  both 
misalignment  parameter  variation  and  aperture- averaged  OPD  determination,  especially  in  the  pm 
regime.  The  need  for  picometer  OPD  tolerances  is  a  relatively  recent  development,  driven  by  ever 
more  demanding  sciece  objectives.  Such  tolerance  requirements  will  likely  become  more 
common,  and  the  lack  of  adequate  analysis  tools  will  correspondingly  be  felt  more  strongly. 


To  develop  a  physical  understanding  of  alignment  sensitivities,  one  would  much  prefer  an 
analytical  rather  than  a  numerical  description  of  the  output  wavefront  as  a  function  of  the 
misalignment  parameters.  Unfortunately,  an  analytical  wavefront  description  of  misaligned 
optical  systems  as  simple  as  a  beam  compressor,  or  even  a  single  focussing  optic,  can  be  too 
complex  to  attempt  by  hand  in  the  kind  of  detail  required  for  sensitivity  studies  (see,  e.g.,  Noecker 
et  al.  1993).  However,  computer  algebra  systems  such  as  Maple  have  advanced  to  such  a  state  of 
capability  and  sophistication  that,  coupled  with  the  processing  power  of  modem  desktop 
computers,  complete  analytical  descriptions  are,  as  we  shall  see,  now  becoming  possible.  This 
paper  describes  AESOP  (An  Extensible  Symbolic  Optics  Package),  an  analytical  ray  tracing 
package  written  in  the  Maple  programming  language. 


AESOP  was  developed  to  support  the  analysis  effort  involved  in  determining  critical  sensitivities 
to  optic  misalignments  in  a  proposed  dual  interferometric  astrometric  telescope,  POINTS 
(Reasenberg  et  al.  1988,  1995a,  1995b,  1996).  POINTS  consists  of  a  pair  of  independent 
Michelson  stellar  interferometers  and  a  laser  metrology  system  that  measures  both  the  critical 
starlight  paths  and  the  angle  between  the  two  interferometer  baselines.  The  nominal  design  has 
baselines  of  2  m,  telescope  apertures  of  35  cm,  and  observes  target  stars  separated  by  roughly  90 
degrees.  One  of  the  distinguishing  features  of  POINTS  is  that  it  employs  holographic  optical 
elements  (HOEs)  to  accomplish  picometer  metrology  over  the  full  aperture  of  the  starlight  optical 
path.  See  the  Reasenberg  et  al.  references  for  a  full  description  of  the  instrument,  its  capabilities, 
and  the  astrophysical,  astrometric,  and  planet-finding  science  that  POINTS  would  significantly 
impact.  Also,  information  can  be  found  on  the  POINTS  web  pages  at 
http://www-cfa.harvard.edu/~reasen/points.html. 


A  key  analysis  problem  regarding  the  POINTS  interferometers  is  the  determination  of  optical  path 
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length  errors  as  a  function  of  various  optical  element  misalignments.  The  path  length  error  budget 
in  a  precision  system  such  as  this  is  only  several  tens  of  picometers.  With  such  a  tight  error 
budget,  it  is  imperative  to  determine  which  perturbations  lead  to  large  path  length  errors.  At  the 
pm  level,  often  we  cannot  trust  our  optical  intuition  in  determining  misalignment  sensitivities.  In 
such  cases,  we  must  rely  on  numerical  analysis  to  an  uncomfortable  degree,  lacking  reliable 
independent  checks  on  the  numerical  results.  AESOP  was  created  in  part  to  fill  this  niche.  In  the 
case  of  POINTS,  a  numerical  program  called  RayTrace  (see  Munson  (1993)  for  a  description)  was 
written  specifically  to  perform  ultra-high  precision,  sub-picometer  OPD  variation  analyses. 
AESOP  was  developed  in  parallel  with  RayTrace.  The  two  analysis  approaches  —  numerical  and 
analytic  —  are  completely  independent  and  therefore  serve  as  excellent  checks  upon  one  another. 


AESOP  traces  an  input  ray  through  a  misaligned  optical  system  and  produces  an  analytic 
description  of  the  output  ray  as  a  function  of  the  system  parameters,  the  misalignment  parameters, 
and  the  input  ray  position  and  direction.  A  crucial  diagnostic  is  the  aperture-averaged  OPD 
variation.  The  physical  principles  involved  are  quite  simple,  since  AESOP  takes  a  classical 
geometric  optics  approach.  At  a  given  reflecting  surface,  all  that  is  required  is  to  calculate  the 
reflected  ray  direction  and  the  accumulated  optical  path  up  to  that  intersection  point.  Similarly,  at 
a  refracting  surface  we  use  Snell's  law  to  calculate  the  refracted  ray  direction.  If  a  holographic 
optical  element  (HOE)  is  encountered,  it  is  a  similarly  simple  process  to  calculate  the  output  ray 
direction  and  the  change  in  optical  phase  across  the  element  (Murison  and  Noecker,  1993).  The 
POINTS  optical  subsystems  involve  all  three  types  of  optical  surfaces.  AESOP  currently  handles 
reflecting  and  holographic  optical  elements.  Refracting  surfaces  will  be  easy  to  add,  due  to  the 
extensible  structure  of  AESOP. 


As  simple  as  the  physics  is,  such  a  ray  tracing  process  is  impossible  to  do  analytically  by  hand 
(especially  aperture- averaged  effects  and  wavefront  analysis).  The  inhibiting  factor  is  the  rapidly 
increasing  (with  each  successive  optical  surface)  complexity  of  the  intermediate  expressions  that 
must  be  algebraically  manipulated.  This  kind  of  repetitive  manipulation  of  unwieldy  objects  is 
precisely  what  computers  can  do  well.  Hence,  a  programmable  computer  algebra  system  like 
Maple  is  well  suited  in  principle  to  analyzing  misaligned  optical  systems,  at  least  simple  ones 
involving  relatively  few  focussing  optical  surfaces. 


In  section  2,  a  few  mathematical  topics  relevant  to  this  kind  of  ray  tracing  in  general  and  to 
AESOP  in  particular  are  briefly  reviewed.  Most  optical  surfaces  are  conicoids,  so  quadric 
surfaces  and  their  normal  vectors  are  introduced  in  section  2.1.  Section  2.2  introduces  the 
concepts  of  optical  rays,  ray  bundles,  and  wavefronts.  Sections  2.3  and  2.4  explain  the  methods 
used  by  AESOP  to  determine  the  intersection  of  a  ray  with  an  optical  surface  and  calculate  the  exit 
ray  direction.  Section  2.5  covers  the  averaging  of  wavefront  error  over  the  aperture  of  a  centrally 
obstructed  optical  system  (common  in  astronomical  systems).  Section  2.6  introduces  the 
expansion  of  a  perturbed  wavefront  in  a  Zemike  series,  the  low-order  components  of  which 
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correspond  to  the  elassical  wavefront  aberrations  such  as  coma,  astigmatism,  etc.  The  Zemike 
series  representation  of  a  wavefront  is  extremely  eonvenient,  instructive,  and  a  helpful  aid  in  the 
analysis  of  perturbed  optical  systems  (as  we  shall  see  in  Appendix  B). 

Assuming  the  baekground  material  in  seetion  2,  seetion  3  eovers  the  AESOP  design  approaeh. 
Optical  systems  and  geometric  ray  tracing  lend  themselves  naturally  to  an  objeet-oriented  design, 
which  AESOP  takes  full  advantage  of.  Indeed,  many  design  simplifications  result,  helping  to 
make  analytical  ray  tracing  not  only  feasible  but  extensible  as  well.  In  practice,  it  has  proven  easy 
to  extend  AESOP  capabilities  as  new  ones  are  needed,  though  one  does  have  to  be  eareful  to 
enforee  the  encapsulation  structure  of  the  AESOP  objeets,  since  the  Maple  language  is  not  itself 
objeet  oriented.  Einally,  section  3.2  gives  an  overview  of  the  AESOP  ray  tracing  process. 
Appendix  A  eontains  a  caleulation  of  the  unit  normal  veetor  for  quadrie  surfaees.  Appendix  B 
ineludes  excerpts  of  a  Maple  session  in  whieh  a  misaligned  beam  eompressor  is  analyzed  by 


AESOP. 


3  2.  Mathematical  Considerations 


3  2.1  Quadric  Surfaces 


3  2.1.1  Surface  Families 


A  general  eonieoid,  or  quadrie  surfaee,  has  the  form 
2 

.  .  P 


s(p)  = 


where  p  is  the  perpendicular  distance  from  the  axis  of  symmetry, /is  the  focal  length  at 


the  vertex,  and  s  is  the  eonie  eonstant.  The  quadrie  surface  types  as  a  funetion  of  e  are: 
e  =  0  paraboloid  of  foeal  length / 
e  =  1  spheroid  of  radius  2 / 

e  <  0  hyperboloid 

0  <  e  <  1  prolate  spheroid 
e  >  1  oblate  spheroid 


The  eccentricity  of  a  prolate  spheroid  is  e  = 


V 


1  -  e 


e 


,  and  that  of  an  oblate  spheroid  is 


e  =  /y - .  Commonly,  the  conic  constant  is  denoted  by  A:  =  e  -  1  so  that  A:  =  0  refers  to  a 

spheroid.  The  quantity  used  here  is  more  convenient  for  our  purposes,  since  the 
paraboloid  is  the  most  frequently  eneountered  focussing  surface  in  astronomical  optics. 
Eor  small  values  of  p  we  have 
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assume ( f >0 ) : 

sag  :=  rho^2  /  (2*f  +  sqrt(4*f^2  -  epsilon*rho^2 ) ) : 

s (rho)  =  series (  sag,  rho,  7  ) ; 


1  1  2  1  £  4  1  6  8 

s(p)  = - p^  + - p^  + - p  +0(p  ) 


The  leading  term  we  recognize  as  the  sagitta  of  a  paraboloid.  All  conicoids  are 
paraboloidal  at  second  order.  Hence,  all  conicoids  exhibit  arbitrarily  good  focussing  of  a 
sufficiently  narrow,  on-axis  input  wavefront.  The  first  aberrational  term  is  spherical 
aberration,  which  enters  in  with  a  linear  term  at  fourth  order  in  radius.  We  see  that 
wavefront  aberrations  are  a  function  of  conicoid  family  type. 


a 


2.1.2  The  Surface  Normal  Vector 


/  2  2 

The  vector  grad  s  is  the  projection  onto  the  xy  plane  (p  =  ^  x  +y  )  of  a  vector  which  is 
normal  to  the  convex  side  of  the  surface  s.  Define  the  surface  function 


g(x,y,  z)  =  z  -  s{x,y).  Then 


A  = 


fa  ^ 

fa  ^ 

,  1 

grad(g) 

_ 

grad(g)  I  I  grad(g) 


is  a  unit  vector  normal  to  the  surface  on  the  concave  side.  Calculation  of  the  surface 
normal  vector  is  the  crucial  step  in  determining  the  output  direction  of  a  ray  incident  on  an 
optical  surface.  As  shown  in  Appendix  A,  we  have  a  relatively  simple  result  for 
conicoids: 


X 

V4/  -  £  P^ 

.  V(l-£)P^  +  4/ 

V(1-£)P^  +  4/ 

V(1-£)P^  +  4/  . 

a 


2.2  Rays,  Ray  Bundles,  and  Wavefronts 


In  the  geometric  optics  regime,  we  may  develop  the  concept  of  a  wavefront  W as  follows. 
Consider  an  infinitely  narrow  beam,  or  ray,  r,  defined  by  an  anchor  point  p,  a  point  in  space 
from  which  the  ray  originates;  a  propagation  direction  v,  conveniently  but  not  necessarily 
represented  as  a  unit  vector;  and  an  optical  path  length  (OPL)  t,  defined  as  the  index  of 
refraction  n  of  the  propagation  medium  integrated  over  a  geometric  distance  L  from  p  along  v: 


t  = 


n( /)  dl.  Hence,  we  may  write  r=p  + tv. 


3 

A  wavefront  may  be  viewed  as  a  surface  W  in  7?  of  constant  optical  phase  propagating 
through  space  (or  through  an  optical  medium).  An  infinitesimally  small  neighborhood  UeW 
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of  each  point peW  propagates  in  a  direction  y(p)  that  is  normal  to  W at p.  We  can  therefore 
associate  a  ray  with  each  point  of  W.  We  define  a  my  bundle  as  the  set  of  rays  belonging  to  W 
.  A  concave  wavefront  (or  portion  of  a  wavefront)  produces  a  converging  ray  bundle,  and  a 
convex  wavefront  produces  a  diverging  ray  bundle. 


In  an  analytical  ray  tracing  procedure,  one  can  consider  a  single  ray  that  is  transformed  by 
passage  through  an  optical  system.  The  resulting  output  ray  then  strikes  a  detector  surface. 
The  position  on  the  detector  of  the  output  ray  intersection  point,  the  output  ray  direction,  and 
the  total  OPL  (the  optical  path  from  the  incident  ray  anchor  point  to  the  detector  surface)  are 
all  functions  of  the  input  ray  anchor  point  and  direction.  Taking  the  input  ray  anchor  point 
position  as  lying  on  an  incident  wavefront  (a  function  we  can  represent  analytically,  usually  a 
plane),  we  can  construct  a  corresponding  output  wavefront  from  the  ray  trace  of  the  input  ray. 


3  2.3  Surface  Intersection  Point 


Given  a  ray  propagating  toward  an  optical  surface,  we  must  find  the  intersection  point  of  the 
ray  with  that  surface.  Define  the  surface  local  coordinate  frame  with  origin  at  the  surface 
vertex  and  Z  axis  along  the  vertex  normal.  An  equation  for  an  input  ray  parameterized  by 
optical  path  t  is 

x{t)  =  r^  +  tv 

where  =  [x^,  y^,  z^]  is  the  ray  anchor  point  and  v  is  the  unit  direction  vector 

V  =  [  sin  \|/  cos  X,  sin  \|/  sin  X,  cos  \|/]  =  [  v  ,  v  ,  v  ] 

X  y  z 

where  (\(/,  X)  are  the  polar  and  azimuthal  angles  with  respect  to  the  +Z  axis  and 
counterclockwise  from  the  +X  axis  in  the  local  coordinate  frame,  respectively.  The  equation 
of  the  surface  —  the  sagitta  —  in  the  local  coordinate  frame  is  of  the  form 

g(x,y,  z)  =  z-s(x.y)  =  0 

We  can  find  the  value  of  t  which  corresponds  to  the  intersection  point  (x,y,z)  by  substituting 

X  =  x^  +  t  sin  \|/  cos  X  =  x^  +  tv 
0  Ox 


y  =  jTq  +  t  sin  \(/  sin  X 
Z  =  Zq  +  t  cos  \|/ 


=  "0  +  ^"z 


into  the  scalar  equation  g(x,y,  z)  and  solving  for  t.  For  a  (perhaps  perturbed)  quadric  surface 
there  will  in  general  be  two  solutions.  AESOP  automatically  chooses  the  correct  solution. 
Then  we  substitute  the  solution  for  t  back  into  the  equations  for  the  ray, 

dt)  =  r^  +  tv 

to  determine  the  x,  y,  and  z  values  of  the  intersection  point. 
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3  2.4  Exit  Ray  Direction 


Upon  encountering  an  optical  surface,  a  ray  must  know  how  to  interact  with  that  surface  and 
choose  an  output  direction.  After  the  intersection  point  is  found,  we  determine  the  unit 
normal  vector  at  that  point,  thus  providing  a  local  reference  for  measuring  input  and  output 
angles.  The  normal  vector  is  easily  found  from  the  gradient  of  the  surface  equation  at  the 
intersection  point. 


a 


2.4.1  Reflection 


For  a  reflecting  surface,  the  incident  and  reflecting  angles  are  equal.  The  reflected  beam 
lies  in  the  plane  defined  by  the  incident  beam  and  the  normal  to  the  surface  at  the 
intersection  point.  Thus,  we  can  write  the  reflected  ray  direction  as  being  equal  to  the 
incident  ray  direction  plus  a  component  along  the  surface  normal  vector  direction. 


V  =v.  +  aA 
r  i 


where  subscripts  r  and  i  correspond  to  the  reflected  and  incident  rays,  respectively,  N  is 
the  normal  to  the  surface  at  the  intersection  point,  and  a  is  a  scale  factor  that  must  be 
determined.  Since  the  incident  and  reflected  angles  are  equal,  we  have 


dot(iV,v  )  dot(iV,v.) 

r  i 

IaIIv  I  |iv||v.| 

r  i 


Combine  this  with 


dot(A,v  )  =  dot(A,  v.l  +  alAp 
r  i 

and  the  further  condition  that  the  magnitude  of  the  reflected  beam  is  equal  to  the 
magnitude  of  the  incident  beam  (certainly  true  if  v.  and  v  are  unit  vectors),  and  we  can 

solve  for  a: 


2dot(A^,  V.) 

i 


a 


After  V  is  determined,  the  result  is  then  transformed  back  to  the  global  reference  frame 
for  propagation  of  the  ray  to  the  next  optical  surface. 

2.4.2  Refraction 


For  a  refracting  surface,  we  must  use  Snell's  law,  which  leads  to  more  complications  than 
the  simple  law  of  reflection.  As  in  the  reflecting  case,  we  may  write 


V  =v.  +  aA 
r  i 


where  the  subscript  r  is  associated  with  the  refracted  ray.  But  now  we  have  the  condition 


n  sin  0  =n  sin  0 
i  i  r  r 
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where  n  and  n  are  the  indices  of  refraction  and  0  and  0  are  the  corresponding  angles  of 
i  r  i  r 

incidence  and  refraction.  Using  dot(A^,  v  )  =  dot(A^,  v.)  +  alA^P  ,  the  condition  I  v  |  =  |v.|, 

o  r  i  r  i 

and  Snell's  law,  we  find  that  the  scale  factor  is 


r 

fn.] 

P  1 

1  1 

i 

■  2  ^ 

V. 

/\  / 1  - 

— 

sin  0  -  cos  0 

1 

n 

i  i 

1  V 

1  n 

y 

where  0  is  determined  from  the  incident  beam  direction  via 

i 


cos  0 .  =  ■ 


dot(A^,  V.) 

i 

livllv  I 

"  i' 


a 


2.4.3  HOE  Diffraction 


a 


2.4. 3.1  Direction  of  the  Diffracted  Ray 


An  adequate  description  of  ray  tracing  across  a  holographic  optical  element  (HOE)  is 
beyond  the  scope  of  this  paper.  However,  it  is  part  of  AESOP's  current  capabilities, 
so  I  present  the  relevant  equations  here,  without  motivation  or  proof.  HOE  ray  tracing 
is  mostly  neglected  in  the  standard  optics  texts,  with  the  exception  of  Welford  (1986). 
Even  in  the  latter,  the  account  is  incomplete,  cursory,  and  potentially  misleading.  The 
reader  is  referred  instead  to  Murison  and  Noecker  (1993)  for  a  complete  and  accurate 
development. 


Define  the  quantity  r=  1  for  transmission  (i.e.,  the  incident  ray  passes  through  the 
HOE),  and  r  =  -1  for  reflection  (for  example  a  diffraction  grating  ruled  on  the  surface 
of  a  mirror).  Then  let  us  define 


S=T  sign(dot(A,  vj) 

where  N  is  now  the  unit  normal  vector  at  the  surface  intersection  point,  and  v .  is  the 

incident  beam  direction  of  propagation  (also  now  required  to  be  a  unit  vector). 
Additionally,  define  the  auxiliary  vector 


m 


X  cross 


N,k 


u  =  cross(A,  V.)  • 

i 


V 


1 


where  m  is  the  diffraction  order  (an  integer),  X  is  the  readout  wavelength  (i.e., 
wavelength  of  the  incident  wavefront),  X  is  the  HOE  construction  wavelength,  and 

k  and  k  are  related  to  the  unit  vectors  directed  from  the  HOE  construction  points 

c .  C- 
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to  the  surface  intersection  point.  If  v, ,  v.,  are  those  unit  vectors,  then  k  =  F,  v,  and 

12  1  1 

k  =  F,  V.,,  where  F ,  =  1  (k=l,2)  if  the  construction  point  C,  is  a  real  focus  and 
1.  z.  k  k 

F  =  -1  if  the  construction  point  C  is  a  virtual  focus.  Then  the  diffracted  ray 

K  K 

direction,  a  unit  vector,  is  given  by 


V  =5 
r 


N-  cross(A^,  u) 


a 


2.43.2  Optical  Path  Correction 


In  geometrical  ray  tracing  of  a  HOE,  a  correction  must  be  added  to  the  optical  path 
upon  traversing  the  HOE  surface.  Again,  refer  to  Murison  and  Noecker  (1993)  for  a 
detailed  development.  The  corrected  optical  path  is 


L  =  L^{p)-M.{p) 

where  ^q(f)  is  the  optical  path  to  the  surface  intersection  point  p  as  calculated  in  the 
normal  geometric  way,  and 


m'k{Xy{p)-Xy{p^)) 

AL(;7)  = - - 

c 

is  the  phase  correction,  p^  is  an  arbitrary  reference  point;  for  a  conicoid,  a  good 

choice  is  the  location  of  the  surface  vertex.  [The  location  of  p^  can  be  arbitrary  since 

it  introduces  a  constant  offset  in  the  optical  path.  We  are  only  interested  in  optical 
path  dijferences  (or  variations).]  The  distance  function  D  is  a  function  of  the  HOE 
construction  points  and  is  given  by  the  expression 


^{p)  =  \V^v^{p)-V^v^{p) 


a 


2.5  OPD  Averaged  over  an  Annular  Aperture 


One  may  define  the  optical  path  difference  (OPD)  as  the  difference  in  total  optical  path 
through  a  system,  starting  from  an  initial  ray  anchor  position  (p,  9),  minus  the  total  optical 
path  of  an  axial  ray  through  the  unperturbed  system  (the  fiducial,  or  chief,  ray).  The  output 
wavefront  is  then  conveniently  represented  by  the  OPD.  Erequently,  we  have  need  of  the 
OPD  averaged  over  the  beam  aperture.  Generally,  there  is  a  central  (usually  circular) 
obscuration,  for  example  the  secondary  mirror  in  a  telescope.  Hence  the  OPD  averaged  over 
an  annular  input  beam  of  inner  and  outer  radii  a  and  b  is 


p  OPD(p,  9)  dd  dp 


OPD 


a  0 


avg 


-a  ) 
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Once  the  OPD  is  determined  by  tracing  a  ray  through  the  system  (and  subtracting  the  fiducial 
ray  optical  path),  the  averaging  integral  is  easy  to  perform.  The  procedure  annular_average() 
is  the  AESOP  function  that  does  this. 


3  2.6  Expansion  of  the  Wavefront  in  a  Zernike  Series 

Zemike  circle  polynomials  are  a  complete  orthogonal  set  over  the  interior  of  the  unit  circle. 
Hence  an  arbitrary  function  W(p,  (|)),  where  p  is  restricted  to  the  range  [0,1],  may  be 
completely  represented  by  an  infinite  series  of  Zemike  polynomials.  We  may  write 


W(p,(^)=  I 
n  =  0 


L  (A  U  (p,(^)  +  5  F  (P,(^)) 
n  n,m  n,m  n,m  n,m 

\m  =  \) 


where  the  values  of  m  are  restricted  to  n-m  =  even,  A  and  B  are  coefficients,  and  U  and  V  are 
given  by 

U  (p,  (|))  =  iJ  I  |(p)  cos(m  (|)) 
n,m  n,\m\ 


V  (p,<^)  =  R 


n,  m 


I  |(P)si 
,\m\ 


sin(w  (|)) 


where  the  radial  polynomials  R  are  given  by 
n  —  m 


R  (p)=  I 
n,m 


h  {n-2k) 

{-\T{n-k)\p 


k\ 


n  +  m 


n  —  m 


See  Murison  (1995)  for  a  discussion,  including  determination  of  the  coefficients  and  an 
example  using  AESOP.  [See  also  Born  and  Wolf  (1980)  and  Zemike  (1934)  for  more 
information  on  Zernike  polynomials.] 


The  Zemike  series  representation  is  useful  for  providing  explicit  expressions  for  the 
well-known  low-order  wavefront  aberrations  such  as  coma,  astigmatism,  defocus,  and  so  on. 
This  turns  out  to  be  an  appealing  way  of  converting  the  often  large  and  inscrutable  AESOP 
wavefront  expressions  into  tidy,  intuitively  understandable  results.  In  general,  the  m=\  terms 
correspond  to  coma,  and  the  m  =  2  terms  correspond  to  astigmatism,  with  n  degrees  of  radial 
"rippliness".  Hence,  the  classical  aberrations  are 


n=l,  m=l 

wavefront  tilt 

n=2,  m=0 

defocus 

n=2,  m=2 

astigmatism 

n=3,  m=l 

coma 

n=4,  m=0 

spherical  aberration 

Another  low-order,  but  non-classical,  aberration  that  is  sometimes  important  is  the  "trefoil" 


Page  10 


term  n=3,  m=3.  Zernike  components  of  the  wavefront  are  illustrated  in  the  example  shown  in 
Appendix  B.  Another  advantage  of  a  Zernike  series  representation  is  that  each  Zernike  term 
affects  the  variance  independently.  Hence,  the  Zernike  polynomials 
minimize  the  wavefront  variance  term  by  term. 


3  3.  AESOP  Design  Considerations 


3  3.1  An  Object  Oriented  Approach 


Geometrical  optics  lends  itself  very  naturally  to  an  object-oriented  approach  when  creating 
computer  programs,  either  numerical  or  algebraic.  AESOP  takes  advantage  of  this  by 
defining  useful  objects  as  Maple  table  structures.  These  Maple  tables  contain,  or  encapsulate, 
all  of  the  information  relevant  to  the  corresponding  objects.  The  Maple  procedures  that 
constitute  AESOP,  and  which  the  user  uses  to  create  a  Maple  procedure  which  can  analyze  an 
optical  system,  manipulate  these  AESOP  objects.  Eollowing  is  a  list  of  the  most  important 
AESOP  objects  with  brief  descriptions. 


3  3.2  AESOP  Objects 


3  3.2.1  Optical  Surface  Data  Structure 


The  optical  system  is  comprised  of  AESOP  optical  elements.  Each  element  type  has  a 
corresponding  procedure  which,  given  certain  information,  creates  the  optical  element 
object.  All  AESOP  optical  elements  share  a  common  table  structure.  The  table  element 
[eqn]  contains  the  equation  describing  the  optical  element's  surface  shape  (usually,  but  not 
necessarily,  a  conicoid)  in  the  surface  local  coordinate  frame.  The  local  frame  origin  is 
located  at  the  surface  vertex,  and  the  positive  Z  axis  is  coincident  with  the  vertex  normal 
vector.  The  [dir]  and  [pos]  elements  contain,  respectively,  a  Maple  vector  and  an  AESOP 
point  which  describe  the  surface  vertex  normal  vector  direction  and  the  vertex  position, 
both  in  the  global  reference  frame.  The  [type]  element  is  the  object  identifier.  Einally,  the 
[coord]  table  element  is  a  point  which  contains  the  (x,  y,  z)  coordinate  labels  that  the  user 
wishes  to  appear  in  the  [eqn]  expression.  An  example  will  make  this  clear: 

read  ' objects. p': 
spheroid{  f,  point { [a, b, c] ) , 


vector { [d [u] , d [v] , d [w] ]) ,  [u,v,w]  ) ; 


tabled 


coord  =  [m,  V,  w] 
type  =  mirror 


flen  =/ 
pos  =[a,b,c] 
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dir  =[d  ,  d  ,  d  ] 

U  V  w 

cc  =  1 

]) 

L  Usually,  one  uses  (x,  y,  z)  for  the  coordinate  labels. 

3  3.2.2  AESOP  Optical  Surface  Objects 


The  current  AESOP  optical  surface  types  are  as  follows.  They  all  have  the  same  Maple 
table  structure  and,  for  the  most  part,  differ  only  in  the  form  of  the  equation  describing  the 
surface.  This  assortment  represents  the  surface  types  needed  in  analyzing  POINTS;  other 
types  are  easy  to  add  as  need  arises. 

3  3.2.2. 1  optical  Jlat 

The  optical_flat  is  a  mirror  with  infinite  focal  length.  It  is  the  simplest  optical 
surface. 


3  3.2.2.2  conicoid 

The  conicoid  is  a  reflecting  surface  with  a  general  conicoid  shape  which  is  a  function 
of  the  conic  constant. 

pos  :=  point ( [a, b, c] ) : 
dir  :=  vector ( [v [x] , v [y] , v [ z ]]) : 
conicoid(  f,  epsilon,  pos,  dir,  [x,y,z]  ) ; 

tabled 

coord =  [x,y,  z] 
type  =  mirror 


eqn  =  z  - 


flen  =/ 

pos  =  [a,  b,  c] 

dir  =  [  V  ,  V  ,  V  ] 
X  y  z 

cc  =  e 

]) 


2  2 
X  +y 


2/+V4/-ex^-e/ 


3  3.2.2. 3  spheroid,  paraboloid 

Because  the  equation  for  a  quadric  surface  simplifies  somewhat  for  the  special  conic 
constant  value  e  =  1,  a  separate  spheroid  surface  is  available.  Similarly,  the 
paraboloid  is  a  conicoid  with  the  special  value  e  =  0. 

3  3.2.2.4  asphere 

The  asphere  object  is  one  whose  conicoidal  surface  is  perturbed  by  a  series  of  radial 
ripples.  AESOP  employs  a  general  asphere  model  of  the  form 
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S(P): 


2/+V4/-ep^ 


I  aJ 


k=  1 


which  is  a  conicoid  plus  a  radial  power  series.  The  AESOP  generating  funetion  for 
this  is  asphereQ.  Here  is  an  illustrative  example: 

clist  :=  [  seq (  A[i],  i=1..4  )  ]: 

asphere (  f,  epsilon,  clist,  pos,  dir,  [x,y,z]  ) ; 

tabled 

coord = [x,y,  z] 
type  =  mirror 


eqn  =  z  - ' 


2  2 
X  +y 


+  A  yj  x^+y^  +A  {x^+y^)+A  (x^+y^) 
2  1  2  3 


3/2 


2f+^4/-ex^-ey 
2 

2  2 
+  A^(xUy^) 

fieri  =f 

pos  =  [a,  b,  c] 

dir  =  [  V  ,  V  ,  V  ] 

X  y  z 

cc  =  e 

]) 

3  3.2.2. 5 pHOE 

The  proeedure  pHOE()  ereates  a  simple  foeussing  HOE  on  a  paraboloidal  mirror  of 
foeal  length  f.  It  is  assumed  that  one  of  the  construetion  points,  say  C^,  is  virtual,  so 

that  a  beam  starting  from  the  other  eonstruction  point,  C^,  will  diffract  to  a  focus  at 
C^.  Eurther  restrictions  are  that  the  diffraetion  order  m  =  1,  and  the  readout 
wavelength  is  equal  to  the  eonstruction  wavelength.  and  must  be  specified  in 
the  surface  local  frame. 

3  3.2.2.6 beam  splitter 

This  is  identieal  to  the  optical_flat  objeet  exeept  for  the  identifieation  tag.  The 
beam_splitter  objeet  exists  solely  for  human  convenience  and  program  readability. 
3  3. 2.2.7 lens 

Currently,  lenses  are  unimplemented.  Deelaration  of  a  lens  will  produce  an  error 
message. 

3  3.2.3  Miscellaneous  Objects 


3  point 

An  AESOP  point  is  identical  in  most  respects  to  a  Maple  vector.  It's  main  purpose  is 
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a 


to  support  the  conceptual  distinction  between  a  direction  vector  and  a  position  (point) 
in  space.  The  AESOP  point  object  is  restricted  to  three  components.  Otherwise,  it  is 
equivalent  to  a  Maple  vector  with  three  elements  and  is  recognized  by  Maple  as  such. 
There  is  a  corresponding  'type/poinT  function  so  that  the  Maple  type  procedures 
(such  as  is(),  hastype(),  etc.)  will  recognize  the  AESOP  point. 
beam 


A  beam  object  is  a  Maple  table  that  has  four  elements.  The  first  two  elements  consist 
of  an  AESOP  point  [pos],  which  contains  the  beam  (or  ray)  anchor  position,  and  a 
Maple  vector  [dir],  which  contains  the  ray  propagation  direction  vector.  Next  is  a 
scalar  element  [path]  for  storing  the  expression  corresponding  to  the  accumulated 
optical  path.  Einally,  a  beam  contains  an  identifier  ['type']  :=  'beam',  which 
procedures  may  query  to  check  that  the  object  is  a  beam.  An  associated  type() 
function  makes  the  Maple  type  procedures  aware  of  beams.  It  is  a  beam  object 
which  serves  as  the  optical  ray  being  propagated  through  an  optical  system. 


a 


3.3  How  AESOP  Does  Ray  Tracing 


An  overview  of  the  ray  tracing  process  using  AESOP  is  as  follows. 


(1)  In  a  Maple  procedure  that  the  user  writes,  the  user  first  defines  the  various  optical  elements 
comprising  the  optical  system.  These  surfaces  are  assembled  into  a  Maple  list  which  AESOP 
routines  will  make  use  of.  Perturbations  (misalignments)  are  applied  in  the  form  of  rotations 
and/or  translations  of  specified  optical  elements.  AESOP  provides  object  rotation  and 
translation  procedures  to  make  this  a  simple  process. 


(2)  The  user  then  defines  the  input  ray,  which  is  subsequently  launched  into  the  optical  system 
by  calling  the  AESOP  procedure  raytrace().  AESOP  then  automatically  traces  the  ray  to  each 
successive  optical  element,  performing  series  expansions  on  the  perturbation  parameter(s)  as 
necessary  and  simplifying  the  cumbersome  expressions  as  as  much  as  possible,  until  finally  an 
output  ray  is  produced  at  the  detector.  Progress  during  this  process  is  communicated  via 
informational  messages  and  key  intermediate  expressions  to  the  monitor  screen.  If  nothing 
else,  there  is  plenty  of  stuff  the  user  can  peruse  while  waiting  for  the  ray  trace  to  finish,  since 
AESOP  is  intentionally  a  bit  chatty. 


(3)  The  OPD  is  then  calculated  from  the  output  ray  expressions,  followed  by  calculation  of  the 
aperture-averaged  OPD. 


(4)  Optionally,  the  Zemike  components  of  the  OPD  are  determined  next,  either  at  the  Maple 
prompt  or  from  within  the  user's  driver  procedure.  The  resulting  Zernike  coefficients  may 
then  be  combined  to  produce  wavefront  aberration  plots.  The  aberrated  wavefronts  are 
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represented  by  3D  Maple  surfaee  plots.  Maple  procedures  are  supplied  for  making  the 
wavefront  plots  in  the  Maple  worksheet. 


An  illustrative  example  of  this  entire  process  for  a  simple  beam  compressor  is  shown  in 
Appendix  B. 


In  practice,  the  essential  step  for  useful  analytical  ray  tracing  is  to  make  series  expansions  at 
each  intersection  of  a  ray  with  an  optical  surface.  (The  original  insight  for  this  trick  is  due  to 
R.D.  Reasenberg.)  This  reduces  the  "equation  bloat"  considerably.  Even  so,  it  is  still  rather 
easy  to  cause  the  intermediate  expressions  to  mushroom  in  size  so  that  they  overwhelm  the 
available  machine  resources.  The  equation  bloat  seems  to  go  as  some  power  of  the  number  of 
focussing  optical  elements  in  a  system.  Flat  surfaces  certainly  contribute  to  increasing 
equation  complexity,  but  at  a  rate  that  pales  in  comparison  to  that  of  focussing  surfaces. 


Since  we  are  interested  in  analyzing  optical  systems  whose  elements  are  slightly  misaligned, 
the  small  parameters  to  perform  the  series  expansion  on  are  naturally  the  misalignment 
perturbations.  Hence,  AESOP  is  not  meant  to  analyze  the  very  interesting  properties  of  ideal, 
perfectly  aligned  optical  systems.  It  requires  at  least  one  misalignment  or  other  perturbation 
parameter. 


For  a  given  optical  system,  a  certain  amount  (sometimes  a  great  amount)  of  tinkering  on  the 
part  of  the  user  is  required  to  hit  upon  the  best  ways  of  simplifying  the  cumbersome 
expressions  so  that  their  size  is  managable.  Great  care  has  been  taken  in  the  types  of 
simplification  taking  place  in  the  AESOP  ray  tracing  routines.  However,  they  are  no  doubt 
optimized  for  the  particular  systems  the  author  has  analyzed  and  will  therefore  perhaps  be  less 
than  optimum  for  other  kinds  of  optical  systems.  Hence,  AESOP  is  nowhere  near  the  "black 
box"  stage,  where  a  user  can  provide  necessary  input,  crank  the  handle,  and  magically  produce 
an  answer  without  caring  overmuch  about  the  internals  of  the  black  box.  Nonetheless, 

AESOP  can  be  quite  useful  and  represents  a  significant  advance  in  capability  for  analyzing 
perturbed  optical  systems  and  performing  misalignment  sensitivity  studies.  It  also  serves  as 
an  invaluable  check  on  numerical  programs  as  well  as  an  essential  aid  to  developing  reliable 
insight  into  the  arcane  and  beautiful  world  of  high-precision  optics. 


3  4.  Availability 

The  AESOP  source  code,  help  files,  examples,  background  papers,  and  other  information  may  be 


found  at  the  web  site  http://aa.usno.navy.mil/AESOP/. 


Page  15 


3  5.  References 


Note:  copies  of  relevant  Smithsonian  Astrophysical  Observatory  (SAO)  Technical  Memoranda 
may  be  obtained  from  the  author.  They  are  also  available  at  http://aa.usno.navy.mil/AESOP/. 


Bom,  M.,  and  Wolf,  E.  (1980).  Principles  of  Optics,  sixth  edition,  Pergammon,  section  9.2. 

Johnston,  K.J.,  Seidelmann,  P.K.,  Germain,  M.E.,  Monet,  D.M.,  Murison,  M.A.,  Urban,  S., 
Davinic,  N.,  Rickard,  E.J.,  Weiler,  K.,  Shao,  M.,  Vaughan,  A.,  Eansen,  J.,  and  Norris,  D.  (1997). 
"Eizeau  Astrometric  Mapping  Explorer  (EAME)",  submitted  to  AstronomicalJournal. 

Eindegren,  E.,  and  Perryman,  M.A.C.  (1996).  "GAIA:  Global  Astrometric  Interferometer  for 
Astrophysics",  ylVron.  and  As  trophy  s.  Supp.  116,  579. 

Eoiseau,  S.,  and  Malbet,  E.  (1996).  "Global  Astrometry  with  OSI",  Astron.  and  Astrophys.  Supp. 
116,  373. 

Eoiseau,  S.,  and  Shaklan,  S.  (1996).  "Optical  Design,  Modelling  and  Tolerancing  of  a  Eizeau 
Interferometer  Dedicated  to  Astrometry",  and  Astrophys.  Supp.  117,  167. 

Murison,  M.A.  (1993).  "Ray  Trace  Analyses  of  Selected  POINTS  Optical  Subsystems",  SAO 
Technical  Memorandum  TM93-08. 

Murison,  M.A.  (1995).  "Expansion  of  Wavefront  Errors  in  an  Infinite  Series  of  Zemike 
Polynomials",  SAO  Technical  Memorandum  TM95-04. 

Murison,  M.A.,  and  Noecker,  M.C.  (1993).  "Ray  Tracing  of  and  Optical  Path  across  a 
Holographic  Optical  Element",  SAO  Technical  Memorandum  TM93-04. 

Noecker,  M.C.,  Murison,  M.A.,  and  Reasenberg,  R.D.  (1993).  "Optic-misalignment  tolerances  for 
the  POINTS  interferometers".  Proceedings  of  the  SPIE  -  The  International  Society  for  Optical 
Engineering,  1947,  218. 

Reasenberg,  R.D.,  Babcock,  R.W.,  Chandler,  J.P.,  Gorenstein,  M.V.,  Huchra,  J.P.,  Pearlman, 

M.R.,  Shapiro,  I.I.,  Taylor,  R.S.,  Bender,  P.,  Buffington,  A.,  Carney,  B.,  Hughes,  J.A.,  Johnston, 
K.I.,  Jones,  B.P.,  and  Matson,  E.E.  (1988).  "Microarcsecond  optical  astrometry  -  An  instmment 
and  its  astrophysical  applications",  Astron.  J.  96,  1731. 

Reasenberg,  R.D.,  Babcock,  R.W.,  Murison,  M.A.,  Noecker,  M.C.,  Phillips,  J.D.,  and  Schumaker, 
B.E.  (1995a).  "POINTS:  The  Instmment  and  its  Mission,"  Proceedings  of  the  SPIE,  Conference 
#2477  on  Spacebome  Interferometry  II  (Orlando,  PE,  17-20  April  1995). 


Page  16 


Reasenberg,  R.D.,  Babcock,  R.W.,  Murison,  M.A.,  Noecker,  M.C.,  Phillips,  J.D.,  Schumaker, 
B.L.,  Ulvestad,  J.S.,  McKinley,  W.,  Zielinski,  R.J.,  and  Lillie,  C.F.  (1995b).  Bull.  American 
Astron.  Soc.,  187,  #71.04. 

R.D.  Reasenberg,  R.W.  Babcock,  M.A.  Murison,  M.C.  Noecker,  J.D.  Phillips,  B.L.  Schumaker, 
J.S.  Ulvestad,  W.  MeKinley,  R.J.  Zielinski,  and  C.F.  Lillie  (1996).  "POINTS:  High  Astrometrie 
Capaeity  at  Modest  Cost  via  Foeused  Design,"  Proceedings  of  the  SPIE,  Conferenee  #2807  on 
Spaee  Teleseopes  and  Instrumentation  IV  (Denver,  CO,  6-7  August  1996),  in  press. 

SIM97.  http://huey.jpl.nasa.gov/sim/ 

Welford,  W.T.  (1986).  Aberrations  of  Optical  Systems,  Adam  Hilger,  Boston. 

Zernike,  F.  (1934).  Physica  1,  689. 


3  6.  Biography 


Marc  A,  Murison  is  an  Astronomer  in  thejAstronomieal  Applieations  Department  6f  the  U.S. 


Naval  Observatory,  in  Washington,  DC.  He  obtained  his  Ph.D.  in  Astronomy  at  the  University  of 
Wisconsin-Madison  in  1988.  He  then  held  a  postdoetoral  position  with  the  Hubble  Spaee 
Teleseope  Wide  Field/Planetary  Camera  team.  Subsequently,  he  was  a  Physieist  at  the 
Smithsonian  Astrophysical  Observatory  where,  among  other  things,  he  was  part  of  the  effort  in 
analyzing  the  optical  systems  of  the  proposed  astrometric  interferometer  satellite,  POINTS.  He 
eurrently  is  part  of  an  effort  by  the  Naval  Observatory  to  build  its  own  astrometrie  satellite.  His 
research  interests  include  the  ehaotie  dynamies  of  the  asteroid  belt,  high-preeision  solar  system 
ephemerides,  and  ultra-high  preeision  analysis  of  optieal  systems. 


a 


Appendix  A:  The  Unit  Normal  Vector  of  a  Quadric  Surface 


The  sagitta  is 

sag; 

2 

P _ 

2/+V4/-ep^ 

Convert  the  sagitta  to  a  funetion  of  (x,y). 

sagxy  :=  subs  (rho=sqrt  (x''2+y''2  ),  sag)  : 

Then  take  the  gradient  and  colleet  on  e. 

map (collect,  grad ( z-sagxy , [x, y , z ] ) ,  [epsilon]): 
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Now  replace  x  and  y  with  p,  where  convenient.  We  get  the  resulting  gradient  vector 

g  :=  map  (  collect,  subs  (x''2+y''2=rho''2 ,  "  )  ,  [x,y,z]  )  ; 


Now  that  we  have  the  gradient  term,  form  the  unit  normal  vector. 

V  :=  array (  [  -dif f ( sagxy , x) ,  -dif f ( sagxy , y ) ,  1  ]  ): 

subs  (  x"'2+y^2=rho^2 ,  evalm(  v  /  mag(g)  )  )  : 

N  :=  map (  collect,  ",  [epsilon, x, y]  ) ; 


2 

P  e 


Vi  1 


2 

P  ye 


2/+V4r-ep 

(f  ( 

1  + 


vv 


(2/+V4r-ep  )  V4r-ep 

2  p^e 
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2/+V4y  -ep  /  2  2,  /  .  2  2 
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where  mag()  is  a  function  which  returns  the  magnitude  of  a  vector.  This  expression  can  be 

2  2  2  2  2 
simplified  considerably.  Let's  get  rid  of  the  x  and  j  terms  and  substitute  Q  for  4y  -  e  p  . 

assume (Q>0 )  : 

map  (  simplify,  N,  {  x^2+y^2=rho^2 ,  4*f ''2-epsilon*rho''2=Q''2  },  [x,y,f]  )  : 

Simplifying  again,  we  get  the  pleasing  result 

map (  simplify,  "  ); 

^  .y  Q 

-  +  'Jp^  +  Q^  Jp^  +  Q^- 

That  is,  we  have 

N  :=  subs (  Q=sqrt (4*f^2-epsilon*rho^2)  ,  "  ): 

N  :=  rootfunc (  N,  collect,  rho  ) ; 


(The  procedure  rootfunc()  applies  an  arbitrary  procedure  to  an  expression  which  contains 
subexpressions  raised  to  fractional  powers.  Now,  the  Maple  procedure  collect()  allows  one  to 
collect  on  a  set  of  variables  and  apply  a  function  to  the  coefficients  of  those  variables,  but  it 
unfortunately  ignores  expressions  raised  to  fractional  powers.  rootfuncQ  also  automatically 
processes  the  components  of  vectors  and  lists.  It  is  a  handy  function  for  easy  manipulation  of 
expressions  under  radicals.) 


a 


Appendix  B:  A  Sample  AESOP  Run 


a 

I 


B.l  OPD  Calculation 
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Following  is  a  Maple  session  in  which  AESOP  is  used  to  calculate  the  aberrated  wavefront  of 
a  misaligned  beam  compressor.  The  primary  mirror  has  focal  length /,  and  the  beam 
compression  ratio  is  denoted  by  C.  The  perturbation  consists  of  a  rotation  of  the  primary 
mirror  about  its  vertex  by  an  angle  0.  The  detector  surface  (a  plane)  is  located  a  distance  d 
below  the  primary  mirror  surface.  The  drive  procedure,  in  which  the  optical  system  is  defined 
and  AESOP  ray  tracing  invoked,  is  called  BeamCompQ.  The  first  argument  of  BeamComp() 
is  the  perturbation  type,  in  this  case  the  primary  mirror  rotation.  The  third  and  fourth 
arguments  are  the  order  of  the  expansions  in  0  and  in  radius. 


read  ' d: /optics/BeamCompressor .p' ; 
#===========================================# 

#  AESOP  # 

#  (An  Extensible  Symbolic  Optics  Processor)  # 

#  version  96.10.02  # 

#===========================================# 

#  Marc  A.  Murison  # 

#  U.S.  Naval  Observatory  # 

#  Astronomical  Applications  Dept.  # 

#  murison@riemann.usno.navy.mil  # 

#  http://aa.usno.navy.mil/AESOP/  # 

#===========================================# 


Use  : 

BeamComp (  pert_typel :: integer , 
epsl : : name, 
pert_type2 : : integer, 
eps2 : : name, 
order_epsl : : integer 
order_eps2 : : integ 
order_r: : integer 

) 

Perturbation  Types 
TRANS_PR1MARY 
ROT_P RIMARY 
TRANS_SECONDARY 
ROT_SECONDARY 
FIELD 
NOTHING 

Example : 

verbosity  : =  1 : 
assume (  Delta  >=  0,  psi  >=  0  )  : 
BeamComp (  TRANS_PR1MARY,  Delta, 


#first  perturbation  type 
#name  of  first  pert,  variable 
tsecond  perturbation  type 
#name  of  second  pert,  variable 
texpansion  order  of  first  pert. 
#expansion  order  of  second  pert. 
#radial  expansion  order 


of  primary  mirror  in  X 
primary  mirror  around  Y  axis 
of  secondary  mirror  in  X 
secondary  mirror  around  Y  axis 
the  input  ray  around  X  axis 


R0T_SEC0NDARY,  psi,  2,  10  ) ; 


translation 
rotation  of 
translation 
rotation  of 
rotation  of 
do  nothing 


NOTE:  you  MUST  assume (epsl  >=  0) :  and  assume (eps2  >=  0) : 
before  calling  BeamComp ( ) 

verbosity  : =  1 : 
assume (  theta  >=  0  ) : 

BeamComp (ROT_PRIMARY, theta, NOTHING,  dummy,  3,12); 

[...much  output  deleted...] 

BeamComp [518] :  Done! 

This  run  took  518  seconds  on  a  100  MHz  Pentium  machine  with  32  MB  of  RAM  available. 
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Peak  memory  usage,  as  reported  by  the  Maple  status  bar,  was  10.4  MB.  The  wavefront  and 
the  aperture- averaged  wavefront,  whieh  are  the  main  results  (and  whieh  are  too  bulky  to 
reproduee  here),  are  stored  in  the  global  variables  OPD  and  OPD_avg. 

3  B.2  Zernike  Series  Decomposition  of  the  Wavefront 

In  this  seetion  we  perform  a  Zernike  series  analysis  of  the  wavefront  (OPD)  just  ealeulated. 
First,  we  set  the  compression  ratio  to  10,  the  distance  to  the  detector  to  20  cm,  and  normalize  p 
so  that  it  spans  the  interval  [0,1]  and  now  R  represents  the  radius  of  the  input  beam.  The  OPD 
simplifies  to 

opd  :=  subs {  rho=R*rho,  C=10,  d=20, 

collect {  OPD,  [theta, cos (phi) , sin  (phi) , rho] ,  simplify  )  ); 


opd 
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1  ( 10940000 +  466841/) 7?  p  1  (51 12000 +  219699/) 7?  p 
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Now  we  perform  the  Zemike  series  expansions  on  the  OPD  to  the  same  order  in  radius  that 
the  ray  trace  calculation  used.  The  procedure  ZernikeSeries()  automatically  calculates  all 
nonzero  Zemike  series  coefficients,  up  to  and  including  the  order  specified.  I  have  clipped  the 
output,  due  to  space  constraints. 

Zemike  Series  (opd,  12,  [theta]  ,  ON,  '  c  '  ,  'S'); 

[...output  deleted...] 

The  Zemike  series  coefficients  are  now  in  a  globally  accessible  table  called  c,  and  the  Zemike 
wavefront  representation  is  stored  in  the  global  variable  S.  Determination  of  the  Zemike 
series  coefficients  is  a  time-consuming  process  (this  particular  example  took  ~11  minutes), 
but  it  is  not  memory  intensive  like  the  ray  tracing  process.  To  show  an  example,  the  classical 
coma  term  is 


ZernikeTerm { 3 , 1 , c) ; 


^  ^  (-2494576896  ^  f  -r  1568023296 


2580480 


+  1568023296  ^  -  57609216000 


V 


818859520  R^  f  +  8962200000  -r  36481536000  19169280000 

61501440000/  + 380813775  268357376o/)  0^  / 

1  R^  {15^  +  1792o/  -r  1344  R^  f  -  320  -  5376  R^  f)  0 


215040 


(3p  -2p)cos((|)) 


We  see  that  it  has  a  first-order  term  in  0  (other  aberrations  begin  at  second  order  in  0).  Hence, 
we  expect  coma  to  be  an  important  aberration  resulting  from  rotational  misalignment  of  the 
primary  mirror.  (In  fact,  for  this  particular  example,  it  turns  out  that  all  of  the  coma  terms 
[3,1],  [5,1],  [7,1],  etc.  are  first  order  in  the  perturbation,  while  the  other  aberrations  are  second 
or  third  order  in  the  perturbation.  Hence,  coma  is  by  far  the  dominant  aberration  after 
wavefront  tilt,  as  we  shall  see  below.) 


We  now  subtract  the  original  OPD  from  the  Zemike  series  representation  to  check  that  the 
two  are  indeed  equivalent. 

factor {  expand { S-opd)  ) ; 

0 

This  is  reassuring!  We  see  by  calculating  the  "cost"  of  each  that  the  Zemike  series 
representation  contains  quite  a  few  more  terms  than  the  heavily  simplified  (in  the  Maple 
sense)  expression  for  OPD. 

cost  { S ) ; 

251  additions  -I-  2433  multiplications  -I-  30  divisions  -I-  17  functions 
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cost (opd) ; 


53  additions  +  599  multiplications  +  28  divisions  +  4  functions 


a 


B.3  Analysis  of  the  Wavefront  Aberrations 


Now  comes  the  fun  part.  We  will  read  in,  among  others,  the  Zemike  plotting  function, 
PlotZernikeWavefrontO-  This  Maple  procedure  produces  a  color  3D  Maple  surface  plot  of 
the  residual  wavefront  after  the  specified  Zemike  terms  have  been  subtracted.  Hence,  it  is  a 
useful  visual  and  quantitative  tool  for  determining  what  are  the  important  aberrations  for  the 
particular  system  and  misalignment  under  study.  Typically,  one  takes  a  look  at  a  plot  of  the 
wavefront,  from  which  the  aberration  with  the  largest  magnitude  is  usually  apparent.  One 
then  subtracts  the  Zemike  term(s)  corresponding  to  that  aberration  and  views  the  resulting 
residual  wavefront,  from  which  the  next  most  important  aberration  is  now  apparent.  This 
process  is  repeated  as  desired,  resulting  in  a  series  of  3D  plots  showing  all  of  the  important 
aberrations.  In  our  example  here,  we  use  a  primary  mirror  focal  length  of  100  cm  and  an  input 
beam  radius  of  10  cm,  which  sets  the  output  beam  radius  at  1  cm  since  we  have  previously 
taken  the  compression  ratio  to  be  10.  For  this  session,  a  perturbation  magnitude  of  0.2  arc 
second  (~1  microradian)  of  primary  mirror  rotation  is  specified.  Changing  the  perturbation 
magnitude  only  changes  the  scale  of  the  aberrations  and  not  their  relative  importance,  as  long 
as  we  remain  in  the  regime  that  is  valid  within  the  expansion  orders  used. 


For  the  first  plot,  we  subtract  the  average  wavefront.  The  resulting  wavefront  residual  is 

verbosity  :=  0: 

PlotZernikeWavefront  (  opd,  theta,  c,  [[0,0]], 

'microns',  evalf (0 . 2*arcsec) ,  100.0,  10.0,  orientation= [-70, 65]  ); 
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Wavefront  Error  (microns) 


10 


(The  fifth  argument  is  a  string  which  sets  the  physical  units  of  the  plot,  in  this  case  microns.) 
The  plot  represents  a  mapping  between  the  input  beam  coordinates  (XY  plane)  and  the  output 
wavefront  (OPD).  Divide  the  horizontal  scale  by  the  compression  factor  (C  =  10  in  this  case) 
to  get  output  beam  coordinates.  We  recognize  that  the  primary  aberration  is,  as  expected, 
wavefront  tilt,  to  the  tune  of  about  1.5  microns  at  the  edge  of  the  beam.  Hence,  let  us 
additionally  remove  the  tilt  term: 

PlotZernikeWavefront (  opd,  theta,  c,  [  [ 0 , 0 ] ,  [  1 , 1 ] ] , 

'microns',  evalf (0 . 2*arcsec) ,  100.0,  10.0  ); 
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Wavefront  Error  (microns) 


The  primary  aberration  is  now  -0.2  microns  of  coma.  All  of  the  coma  terms  are  first  order  in 
the  rotation  angle  0,  and  therefore  dominant.  Let  us  then  additionally  subtract  all  orders  of 
coma: 

PlotZernikeWavef ront (  opd,  theta,  c, 

'pm',  evalf ( 0 . 2*arcsec) ,  100.0,  10.0  ); 
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Wavefront  Error  (pm) 


80 
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OPD 

40 
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.  0 


radius  (cm)  -5  1 0 


The  residuals  now  consist  of  a  -100  picomeX&r  mix  of  defocus  and  astigmatism,  which  both 
are  quadratic  in  radius.  By  removing  all  of  the  first  order  (in  0)  terms,  we  see  that  the  second 
order  terms  result  in  aberrations  that  are  over  2,000  times  smaller.  Similar  to  the  series  of 
coma  terms,  whose  higher-order  (in  radius)  components  were  of  similar  magnitude  to  the 
classical  third  order  coma  term  [3,1],  the  defocus  and  astigmatism  terms  of  various  radial 
orders  contribute  comparable  amounts  to  the  aberrations  as  the  classical  components  [2,0]  and 
[2,2],  respectively.  Hence,  let  us  additionally  remove  all  defocus  and  astigmatism  terms  to  get 

PlotZernikeWavef ront {  opd,  theta,  c, 

[[0,0],  [1,1],  [3,1],  [5,1],  [7,1],  [9,1],  [11,1],  [2,0],  [2, 2],  [4,0], 

[4,2],  [6,0],  [6,2],  [8,0],  [8,2],  [10,0],  [10,2],  [12,0],  [12,2]], 

' fm' ,  evalf (0 . 2*arcsec) ,  100.0,  10.0  ); 


Page  26 


Wavefront  Error  (fm) 


0.2 

O.V 

OPD  0 

-0.1 

-0. 

10 


radius  (cm) 


-1010 


Now  we  recognize  exceedingly  small  0.2  fm  residuals  due  to  trefoil  ([3,3],  [5,3],  etc.) 
aberrations.  These  aberrations  happen  in  this  system  to  be  third  order  in  the  perturbation 
angle  0.  For  example, 

ZernikeTerm ( 3 , 3, c) ; 

— 0^  (49018305  ^  f+  1 148700000  -  105455520  R^  f"  -  2453760000  R^  ? 

860160 

-t  200878720  R^ ^  +  4632320000  R^  f  -  315227136  R^/  -  7182336000  R^/ 

+  337559040/  +  7526400000/)p^cos(3  4))  / 


a 


B.4  Conclusion 


We  have  firmly  (and  simply!)  established  that  this  particular  wavefront's  dominant  aberrations 
are  tilt  and  coma,  followed  at  a  much  lower  level  by  defocus  and  astigmatism.  Even  more 
important,  we  have  the  dependence  of  each  aberration  type  (as  well  as  of  the 
aperture- averaged  wavefront)  on  the  perturbation  parameter  0,  as  well  as  the  optical  system 
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parameters,  allowing  us  to  determine  the  sensitivity  of  these  aberrations  to  (in  this  case) 
rotation  of  the  primary  mirror.  Hence,  given  an  error  budget,  one  can  easily  determine 
tolerances  for  the  perturbation  magnitudes. 

ave  OPD,  OPD_avg,  opd,  S,  c,  ' d : /optics/AESOP/MapleTech/AESOP_article . eqs ' ; 
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